hide
# Set ggplot theme for figures
theme_set(theme_bw() +
theme(
# Titles
plot.title = element_text(hjust = 0, face = "bold", size = 14),
plot.subtitle = element_text(size = 9),
plot.caption = element_text(size = 7.2, margin = margin(t = 20), color = "gray40"),
legend.title = element_text(size = 9),
legend.text = element_text(size = 7.5),
# Axes
axis.line.y = element_line(),
axis.ticks.y = element_line(),
axis.line.x = element_line(),
axis.ticks.x = element_line(),
axis.text = element_text(size = 11),
# Facets
strip.text = element_text(color = "white",
face = "bold",
size = 11),
strip.background = element_rect(
color = "white",
fill = "#00736D",
size = 1,
linetype="solid"))
)
# Building a GMRI theme based on Wall street Journal and NYTimes theme
# base settings from {ggthemes}
theme_gmri <- function(base_size = 10,
bg_color = "lightblue",
base_family = "sans",
title_family = "sans",
facet_color = "teal") {
# Color from gmRi palette, sets background color
#colorhex <- gmRi::gmri_cols()[bg_color]
facet_hex <- gmri_cols()[facet_color]
# Set up theme
theme_foundation(
base_size = base_size,
base_family = base_family) +
theme(
# Major Elements
line = element_line(linetype = 1, colour = "black"),
rect = element_rect(fill = "transparent",
linetype = 0,
colour = NA),
text = element_text(colour = "black"),
title = element_text(family = title_family, size = 12),
# Axis elements
axis.text.x = element_text(colour = NULL),
axis.text.y = element_text(colour = NULL),
axis.ticks = element_line(colour = NULL),
axis.ticks.y = element_blank(),
axis.ticks.x = element_line(colour = NULL),
axis.line = element_line(),
axis.line.y = element_blank(),
axis.text = element_text(size = 11),
# Legend Elements
legend.background = element_rect(),
legend.position = "top",
legend.direction = "horizontal",
legend.box = "vertical",
legend.title = element_text(size = 9),
legend.text = element_text(size = 7.5),
# Panel/Grid Setup
panel.grid = element_line(colour = NULL, linetype = 3, color = "gray70"),
panel.grid.major = element_line(colour = "black"),
panel.grid.major.x = element_blank(),
panel.grid.minor = element_blank(),
# Title and Caption Details
plot.title = element_text(hjust = 0, face = "bold", size = 14),
plot.subtitle = element_text(size = 9),
plot.caption = element_text(size = 7.2, margin = margin(t = 20)),
plot.margin = unit(c(1, 1, 1, 1), "lines"),
# Facet Details
strip.text = element_text(color = "white", face = "bold", size = 11),
strip.background = element_rect(
color = "white",
fill = facet_hex,
size = 1,
linetype="solid"))
}
# Set theme up for maps
map_theme <- list(
theme(
panel.border = element_rect(color = "black", fill = NA),
plot.background = element_rect(color = "transparent", fill = "transparent"),
line = element_blank(),
axis.title.x = element_blank(), # turn off titles
axis.title.y = element_blank(),
#legend.position = "bottom",
legend.title.align = 0.5))
2021 Sea Surface Temperature
An Exceptional Year
2021 was an exceptionally warm year for the Gulf of Maine. It had the warmest fall on record for the area, and the second warmest summer. Much of the year the region experienced what would be considered “heatwave conditions”.
hide
# Load the timeseries
region_timeseries <- read_csv(timeseries_path, col_types = cols(), guess_max = 1e6)
# Clean up the data - add labels
region_timeseries <- region_timeseries %>%
distinct(time, .keep_all = T) %>%
mutate(time = as.Date(time),
yr = year(time),
month_num = month(time),
month = month(time, label = T, abbr = T),
week = lubridate::week(time),
doy = yday(time),
season = metR::season(month_num, lang = "en"),
#Set up correct year for winters, they carry across into next year
season_eng = case_when(
season == "SON" ~ "Fall",
season == "DJF" ~ "Winter",
season == "MAM" ~ "Spring",
season == "JJA" ~ "Summer"),
season_yr = ifelse( (season_eng == "Winter" & month_num %in% c(1,2)), yr - 1, yr)
)
# Get heatwave statuses for each day:
# Uses area weighted sst by default
region_hw <- pull_heatwave_events(
temperature_timeseries = region_timeseries,
threshold = 90,
clim_ref_period = c("1982-01-01", "2011-12-31"))
# Format dates and seasons and heatwave outcomes
region_hw <- region_hw %>%
mutate(time = as.Date(time),
yr = year(time),
month_num = month(time),
month = month(time, label = T, abbr = T),
week = lubridate::week(time),
doy = yday(time),
season = metR::season(month_num, lang = "en"),
season_eng = case_when(
season == "SON" ~ "Fall",
season == "DJF" ~ "Winter",
season == "MAM" ~ "Spring",
season == "JJA" ~ "Summer"),
#Set up correct year for winters, they carry across into next year
season_yr = ifelse(
(season_eng == "Winter" & month_num %in% c(1,2)), yr - 1, yr)
)
# Set up heatwave status outcomes for fancy table features:
region_hw <- region_hw %>%
mutate(hw_outcomes = case_when(
mhw_event == TRUE ~ 1,
mcs_event == TRUE ~ 0,
TRUE ~ 0.5))
Departure from Historic Conditions
Thanks to some long-running temperature monitoring programs we are able to track monthly sea surface temperatures as far back as 1850. The figure below shows how the region’s temperatures have fluctuated as measured by the Extended Reconstructed Sea Surface Temperature (ERSST) record.
Once we hit 1982 the availability of a higher resolution temperature resource becomes available in the form of the Optimum Interpolation Sea Surface Temperature (OISST) record. The ERSST data is of a coarser resolution (monthly measurements at a 0.5 x 0.5 degree grid) and does not capture the warmer inshore dynamics (a bias to show colder temperature), but is within half of a degree of the more modern measurements.
hide
# Pathto ERSST on box
ersst_path <- cs_path("okn", "ersst")
# Load temps and anoms (1982-2011 climatology)
ersst_temp <- stack(str_c(ersst_path, "ERSSTv5_merged.nc"))
ersst_anom <- stack(str_c(ersst_path, "ERSSTv5_anom.nc"))
# Crop GoM
# function to mask
mask_nc <- function(ras_obj, mask_shape){
# Get stats from ranks
m1 <- raster::mask(rotate(ras_obj), mask_shape)
m1 <- raster::crop(m1, mask_shape)
return(m1)}
# Load the shape, mask ERSST with it
gom_shape <- read_sf(poly_path)
gom_ersst <- mask_nc(ras_obj = ersst_temp, mask_shape = gom_shape)
gom_ersst_anom <- mask_nc(ras_obj = ersst_anom, mask_shape = gom_shape)
# Clean them up - temp
ersst_temp_tl <- as.data.frame(cellStats(gom_ersst, mean)) %>%
rownames_to_column(var = "Date") %>%
setNames(c("Date", "sst")) %>%
mutate(Date = str_sub(Date, 2, -1),
Date = str_replace_all(Date, "[.]", "-"),
Date = as.Date(Date),
yr = str_sub(Date, 1, 4))
# anomalies
ersst_anom_tl <- as.data.frame(cellStats(gom_ersst_anom, mean)) %>%
rownames_to_column(var = "Date") %>%
setNames(c("Date", "sst_anom")) %>%
mutate(Date = str_sub(Date, 2, -1),
Date = str_replace_all(Date, "[.]", "-"),
Date = as.Date(Date),
yr = str_sub(Date, 1, 4))
# Join them together
ersst_tl <- left_join(ersst_temp_tl, ersst_anom_tl) %>%
mutate(yr = as.numeric(yr))
# Get Yearly Average
ersst_yr <- ersst_tl %>%
group_by(yr) %>%
summarise(sst = mean(sst),
sst_anom = mean(sst_anom),
.groups = "drop") %>%
mutate(yr = as.numeric(yr)) %>%
filter(yr <= 2019)
# Same for OISST
oisst_yr <- region_hw %>%
filter(between(yr, 1982, 2021)) %>%
group_by(yr) %>%
summarise(sst = mean(sst), .groups = "drop")
# Make Splines
ersst_smooth <- as.data.frame(spline(ersst_yr$yr, ersst_yr$sst))
oisst_smooth <- as.data.frame(spline(oisst_yr$yr, oisst_yr$sst))
# Plot them
ggplot() +
geom_line(data = ersst_smooth, aes(x, y, color = "ERSSTv5"),
group = 1, alpha = 0.4, linetype = 1) +
geom_point(data = ersst_yr, aes(yr, sst, color = "ERSSTv5"),
size = 1) +
geom_line(data = oisst_smooth, aes(x, y, color = "OISSTv2"),
group = 1, alpha = 0.4, linetype = 1) +
geom_point(data = oisst_yr, aes(yr, sst, color = "OISSTv2"),
size = 1) +
scale_color_gmri() +
theme_gmri() +
labs(x = "", y = "Sea Surface Temperature",
color = "Temperature Data Record",
title = "Long-Term Temperature Record for the Gulf of Maine",
subtitle = "Current temperatures above 150-year highs") +
scale_y_continuous(sec.axis = sec_axis(trans = ~as_fahrenheit(., data_type = "anomalies"),
labels = number_format(suffix = " \u00b0F")),
labels = number_format(suffix = " \u00b0C")) +
guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) +
facet_zoom(xlim = c(1982, 2021)) +
theme(legend.position = "bottom")
Relationship to Climate Patterns
To better get a sense of what kind of year 2021 was it is helpful to relate it to what we would expect from similar years.
ENSO Index
One large feature in global weather patterns is the El Niño Southern Oscillation (ENSO). While ENSO and the southern oscillation index are derived from temperature/pressure differences in the Pacific Ocean, they influence atmospheric behavior coming across the United States which has the potential to influence SST patterns.
hide
# Source: https://www.ncdc.noaa.gov/teleconnections/enso/soi
# Does not work:
# library(rsoi)
# clim_idx <- download_soi()
# Manually copied from link above
soi_std <- tibble::tribble(#####
~"YEAR", ~"JAN", ~"FEB", ~"MAR", ~"APR", ~"MAY", ~"JUN", ~"JUL", ~"AUG", ~"SEP", ~"OCT", ~"NOV", ~"DEC",
1951, 1.5, 0.9, -0.1, -0.3, -0.7, 0.2, -1.0, -0.2, -1.1, -1.0, -0.8, -0.7,
1952, -0.9, -0.6, 0.5, -0.2, 0.8, 0.7, 0.5, 0.1, -0.2, 0.4, 0.0, -1.2,
1953, 0.3, -0.5, -0.2, 0.2, -1.7, 0.1, -0.0, -1.2, -1.2, 0.1, -0.3, -0.5,
1954, 0.7, -0.3, 0.3, 0.6, 0.5, 0.1, 0.4, 1.1, 0.2, 0.3, 0.1, 1.4,
1955, -0.5, 1.9, 0.6, -0.1, 1.0, 1.3, 1.6, 1.5, 1.3, 1.5, 1.2, 1.0,
1956, 1.3, 1.6, 1.3, 0.9, 1.4, 1.1, 1.1, 1.2, 0.1, 1.8, 0.2, 1.1,
1957, 0.6, -0.1, 0.2, 0.2, -0.7, 0.2, 0.2, -0.5, -0.9, 0.0, -1.0, -0.3,
1958, -1.9, -0.5, 0.3, 0.4, -0.5, 0.3, 0.4, 0.9, -0.3, 0.1, -0.4, -0.6,
1959, -0.9, -1.4, 1.3, 0.4, 0.5, -0.1, -0.3, -0.1, 0.0, 0.5, 0.9, 0.9,
1960, 0.1, 0.1, 1.0, 0.8, 0.5, 0.1, 0.5, 0.8, 0.7, 0.1, 0.5, 0.8,
1961, -0.3, 0.9, -1.8, 0.8, 0.3, 0.1, 0.2, 0.2, 0.1, -0.3, 0.5, 1.5,
1962, 2.0, -0.3, 0.1, 0.2, 1.1, 0.7, 0.1, 0.6, 0.4, 1.0, 0.3, 0.2,
1963, 1.0, 0.6, 1.1, 0.8, 0.4, -0.5, -0.1, 0.0, -0.6, -1.2, -0.8, -1.2,
1964, -0.4, 0.0, 1.1, 1.1, 0.2, 0.8, 0.6, 1.5, 1.3, 1.3, 0.2, -0.3,
1965, -0.4, 0.4, 0.8, -0.5, 0.2, -0.6, -1.8, -0.7, -1.3, -0.9, -1.5, 0.2,
1966, -1.3, -0.2, -0.9, -0.2, -0.4, 0.3, 0.1, 0.6, -0.2, -0.1, -0.0, -0.3,
1967, 1.7, 1.7, 1.2, -0.0, -0.0, 0.6, 0.2, 0.7, 0.5, 0.1, -0.4, -0.6,
1968, 0.5, 1.3, 0.1, 0.0, 1.2, 1.1, 0.7, 0.3, -0.3, -0.1, -0.3, 0.2,
1969, -1.5, -0.5, 0.4, -0.4, -0.2, 0.2, -0.5, -0.1, -1.0, -0.9, -0.1, 0.4,
1970, -1.1, -1.0, 0.6, -0.1, 0.4, 1.0, -0.4, 0.6, 1.2, 1.0, 1.6, 1.9,
1971, 0.4, 2.0, 2.3, 1.7, 0.9, 0.4, 0.2, 1.5, 1.4, 1.7, 0.5, 0.3,
1972, 0.5, 1.1, 0.6, -0.1, -1.6, -0.5, -1.4, -0.5, -1.4, -0.9, -0.3, -1.3,
1973, -0.3, -1.4, 0.7, 0.1, 0.4, 1.1, 0.6, 1.3, 1.2, 0.8, 2.6, 1.8,
1974, 2.4, 2.1, 2.4, 0.9, 1.0, 0.4, 1.1, 0.8, 1.1, 0.9, -0.1, 0.2,
1975, -0.5, 0.8, 1.6, 1.2, 0.6, 1.3, 1.9, 2.0, 2.1, 1.7, 1.2, 2.1,
1976, 1.4, 1.7, 1.7, 0.3, 0.4, 0.3, -0.9, -0.8, -1.1, 0.4, 0.7, -0.3,
1977, -0.4, 1.2, -0.5, -0.4, -0.5, -0.9, -1.1, -0.8, -0.8, -1.0, -1.3, -1.1,
1978, -0.3, -2.7, -0.2, -0.3, 1.4, 0.7, 0.6, 0.4, 0.1, -0.4, -0.0, -0.1,
1979, -0.4, 1.0, 0.1, -0.1, 0.5, 0.6, 1.3, -0.2, 0.1, -0.1, -0.4, -0.7,
1980, 0.4, 0.3, -0.4, -0.6, -0.0, -0.0, -0.0, 0.4, -0.5, 0.0, -0.3, -0.1,
1981, 0.4, -0.2, -1.3, -0.1, 0.8, 1.2, 0.8, 0.7, 0.3, -0.4, 0.2, 0.5,
1982, 1.2, 0.3, 0.6, 0.1, -0.3, -1.0, -1.5, -1.7, -1.7, -1.7, -2.6, -2.2,
1983, -3.5, -3.6, -2.4, -0.9, 0.6, 0.0, -0.6, 0.1, 0.9, 0.4, -0.1, 0.0,
1984, 0.2, 0.9, -0.2, 0.3, 0.2, -0.3, 0.2, 0.4, 0.1, -0.3, 0.3, -0.1,
1985, -0.3, 1.2, 0.8, 1.2, 0.4, -0.4, -0.1, 1.0, 0.0, -0.4, -0.2, 0.2,
1986, 1.0, -1.0, 0.5, 0.3, -0.2, 1.0, 0.3, -0.4, -0.5, 0.6, -1.2, -1.4,
1987, -0.7, -1.2, -1.3, -1.4, -1.3, -1.1, -1.4, -0.9, -1.0, -0.4, -0.0, -0.5,
1988, -0.1, -0.4, 0.6, 0.1, 0.9, 0.1, 1.0, 1.5, 1.8, 1.4, 1.7, 1.2,
1989, 1.5, 1.2, 1.1, 1.6, 1.2, 0.7, 0.9, -0.3, 0.5, 0.8, -0.2, -0.5,
1990, -0.1, -1.8, -0.4, 0.2, 1.2, 0.3, 0.5, -0.2, -0.7, 0.3, -0.5, -0.2,
1991, 0.6, 0.3, -0.7, -0.6, -1.0, -0.1, 0.0, -0.4, -1.5, -1.0, -0.7, -1.8,
1992, -2.9, -0.9, -2.0, -1.0, 0.3, -0.6, -0.6, 0.4, 0.1, -1.4, -0.7, -0.6,
1993, -0.9, -0.7, -0.5, -1.2, -0.3, -0.8, -0.8, -0.9, -0.7, -1.1, -0.1, 0.2,
1994, -0.1, 0.3, -0.7, -1.3, -0.7, -0.4, -1.3, -1.2, -1.6, -1.1, -0.6, -1.2,
1995, -0.4, -0.1, 0.8, -0.7, -0.4, 0.1, 0.4, 0.3, 0.3, 0.0, 0.0, -0.5,
1996, 1.0, 0.3, 1.1, 0.8, 0.3, 1.2, 0.7, 0.7, 0.6, 0.6, -0.1, 0.9,
1997, 0.5, 1.7, -0.4, -0.6, -1.3, -1.4, -0.8, -1.4, -1.4, -1.5, -1.2, -1.0,
1998, -2.7, -2.0, -2.4, -1.4, 0.3, 1.0, 1.2, 1.2, 1.0, 1.1, 1.0, 1.4,
1999, 1.8, 1.0, 1.3, 1.4, 0.2, 0.3, 0.5, 0.4, -0.1, 1.0, 1.0, 1.4,
2000, 0.7, 1.7, 1.3, 1.2, 0.4, -0.2, -0.2, 0.7, 0.9, 1.1, 1.8, 0.8,
2001, 1.0, 1.7, 0.9, 0.2, -0.5, 0.3, -0.2, -0.4, 0.2, -0.0, 0.7, -0.8,
2002, 0.4, 1.1, -0.2, -0.1, -0.8, -0.2, -0.5, -1.0, -0.6, -0.4, -0.5, -1.1,
2003, -0.2, -0.7, -0.3, -0.1, -0.3, -0.6, 0.3, 0.1, -0.1, 0.0, -0.3, 1.1,
2004, -1.3, 1.2, 0.4, -0.9, 1.0, -0.8, -0.5, -0.3, -0.3, -0.1, -0.7, -0.8,
2005, 0.3, -3.1, 0.3, -0.6, -0.8, 0.4, 0.2, -0.3, 0.4, 1.2, -0.2, -0.0,
2006, 1.7, 0.1, 1.8, 1.1, -0.5, -0.2, -0.6, -1.0, -0.6, -1.3, 0.1, -0.3,
2007, -0.8, -0.1, 0.2, -0.1, -0.1, 0.5, -0.3, 0.4, 0.2, 0.7, 0.9, 1.7,
2008, 1.8, 2.6, 1.4, 0.7, -0.1, 0.6, 0.3, 1.0, 1.2, 1.3, 1.3, 1.4,
2009, 1.1, 1.9, 0.4, 0.8, -0.1, 0.1, 0.2, -0.2, 0.3, -1.2, -0.6, -0.7,
2010, -1.1, -1.5, -0.7, 1.2, 0.9, 0.4, 1.8, 1.8, 2.2, 1.7, 1.3, 2.9,
2011, 2.3, 2.7, 2.5, 1.9, 0.4, 0.2, 1.0, 0.4, 1.0, 0.8, 1.1, 2.5,
2012, 1.1, 0.5, 0.7, -0.3, 0.0, -0.4, -0.0, -0.2, 0.2, 0.3, 0.3, -0.6,
2013, -0.1, -0.2, 1.5, 0.2, 0.8, 1.2, 0.8, 0.2, 0.3, -0.1, 0.7, 0.1,
2014, 1.4, 0.1, -0.9, 0.8, 0.5, 0.2, -0.2, -0.7, -0.7, -0.6, -0.9, -0.6,
2015, -0.8, 0.2, -0.7, -0.0, -0.7, -0.6, -1.1, -1.4, -1.6, -1.7, -0.5, -0.6,
2016, -2.2, -2.0, -0.1, -1.2, 0.4, 0.6, 0.4, 0.7, 1.2, -0.3, -0.1, 0.3,
2017, 0.2, -0.1, 0.9, -0.2, 0.3, -0.4, 0.8, 0.5, 0.6, 0.9, 0.9, -0.1,
2018, 1.1, -0.5, 1.5, 0.5, 0.4, -0.1, 0.2, -0.3, -0.9, 0.4, -0.1, 1.0,
2019, -0.0, -1.4, -0.3, 0.1, -0.4, -0.5, -0.4, -0.1, -1.2, -0.4, -0.8, -0.6,
2020, 0.2, -0.1, -0.1, 0.2, 0.4, -0.4, 0.4, 1.1, 0.9, 0.5, 0.7, 1.8,
2021, 1.9, 1.5, 0.4, 0.3, 0.5, 0.4, 1.4, 0.6, 0.8, 0.7, 1.0, 1.5
)
#####
# Pick a date centered in each month so you can plot them as one timeline
month_centers <- data.frame(Month = toupper(month.abb),
Month_num = str_pad(seq(1,12,1), width = 2, pad = "0", side = "left"))
# Reformat
soi_long <- soi_std %>%
pivot_longer(names_to = "Month", values_to = "soi_z", cols = JAN:DEC) %>%
left_join(month_centers, by = "Month") %>%
mutate(Date = as.Date(str_c(YEAR, "-", Month_num, "-15")),
ribbon_fill = ifelse(soi_z > 0,"above", "below"),
ribbon_above = ifelse(soi_z > 0, soi_z, NA),
above_ymin = ifelse(soi_z > 0, 0, NA),
ribbon_below = ifelse(soi_z < 0, soi_z, NA),
below_ymin = ifelse(soi_z < 0, 0, NA)
)
# Set color scale and break points
cutpoints <- rev(seq(-4, 4, 1))
col_rdbu <- rev(RColorBrewer::brewer.pal(n = length(cutpoints), name = "RdBu"))
# fix GeomRibbon
GeomRibbon$handle_na <- function(data, params) { data }
# Make figure
ggplot(filter(soi_long, YEAR >= 1982), aes(x = Date, y = soi_z)) +
geom_ribbon(aes(ymin = above_ymin, ymax = ribbon_above, fill = "La Niña"), alpha = 0.8) +
geom_ribbon(aes(ymin = below_ymin, ymax = ribbon_below, fill = "El Niño"), alpha = 0.8) +
scale_fill_manual(values = c("La Niña" = col_rdbu[1], "El Niño" = col_rdbu[8])) +
geom_line(size = 0.25, color = "gray20") +
labs(x = "Date", y = "SOI", title = "ENSO Index",
caption = "Years overlapping OISST data displayed", fill = "") +
theme_gmri() +
theme(legend.position = "bottom")
When we look at average temperatures for the Gulf of Maine, and how that relates to the ENSO index, we found there was NO relationship. Meaning that a current month’s El niño/ La Nińa status did not correlate with warmer/hotter months.
North Atlantic Oscillation
A similar climate index to ENSO, with more relevance to our area is the North Atlantic Oscillation (NAO). The NAO index is based on the surface sea-level pressure difference between the Subtropical (Azores) High and the Subpolar Low.
Strong positive phases of the NAO tend to be associated with above-normal temperatures in the eastern United States and across northern Europe and below-normal temperatures in Greenland and oftentimes across southern Europe and the Middle East.
hide
# NAO Data Source:
# https://www.ncdc.noaa.gov/teleconnections/nao/
nao_idx <- tribble( #####
~"YEAR", ~"Jan", ~"Feb", ~"Mar", ~"Apr", ~"May", ~"Jun", ~"Jul", ~"Aug", ~"Sep", ~"Oct", ~"Nov", ~"Dec",
1950, 0.92, 0.40, -0.36, 0.73, -0.59, -0.06, -1.26, -0.05, 0.25, 0.85, -1.26, -1.02,
1951, 0.08, 0.70, -1.02, -0.22, -0.59, -1.64, 1.37, -0.22, -1.36, 1.87, -0.39, 1.32,
1952, 0.93, -0.83, -1.49, 1.01, -1.12, -0.40, -0.09, -0.28, -0.54, -0.73, -1.13, -0.43,
1953, 0.33, -0.49, -0.04, -1.67, -0.66, 1.09, 0.40, -0.71, -0.35, 1.32, 1.04, -0.47,
1954, 0.37, 0.74, -0.83, 1.34, -0.09, -0.25, -0.60, -1.90, -0.44, 0.60, 0.40, 0.69,
1955, -1.84, -1.12, -0.53, -0.42, -0.34, -1.10, 1.76, 1.07, 0.32, -1.47, -1.29, 0.17,
1956, -0.22, -1.12, -0.05, -1.06, 2.21, 0.10, -0.75, -1.37, 0.24, 0.88, 0.51, 0.10,
1957, 1.05, 0.11, -1.26, 0.49, -0.79, -0.72, -1.19, -0.55, -1.66, 1.32, 0.73, 0.12,
1958, -0.54, -1.06, -1.96, 0.37, -0.24, -1.38, -1.73, -1.56, -0.07, 0.16, 1.64, -0.70,
1959, -0.87, 0.68, -0.15, 0.36, 0.39, 0.40, 0.74, 0.06, 0.88, 0.89, 0.41, 0.44,
1960, -1.29, -1.89, -0.50, 1.36, 0.45, -0.21, 0.35, -1.40, 0.39, -1.73, -0.51, 0.06,
1961, 0.41, 0.45, 0.55, -1.55, -0.36, 0.86, -0.39, 0.90, 1.24, 0.51, -0.62, -1.48,
1962, 0.61, 0.55, -2.47, 0.99, -0.10, 0.16, -2.47, 0.14, -0.37, 0.41, -0.23, -1.32,
1963, -2.12, -0.96, -0.43, -1.35, 2.16, -0.43, -0.77, -0.64, 1.79, 0.94, -1.27, -1.92,
1964, -0.95, -1.43, -1.20, 0.36, 0.52, 1.29, 1.90, -1.77, 0.20, 0.74, -0.01, -0.15,
1965, -0.12, -1.55, -1.51, 0.72, -0.62, 0.29, 0.32, 0.45, 0.37, 0.38, -1.66, 1.37,
1966, -1.74, -1.39, 0.56, -0.75, 0.22, 1.05, 0.32, -1.76, -0.45, -0.68, -0.04, 0.72,
1967, -0.89, 0.19, 1.51, 0.18, -0.99, 1.40, 0.41, 1.44, 0.93, 0.07, 0.60, -0.45,
1968, 0.13, -1.29, 0.40, -1.08, -1.76, 0.33, -0.80, -0.66, -1.92, -2.30, -0.93, -1.40,
1969, -0.83, -1.55, -1.56, 1.53, 0.55, 0.55, 0.57, -1.45, 2.07, 0.66, -0.96, -0.28,
1970, -1.50, 0.64, -0.96, -1.30, 1.14, 1.55, 0.10, 0.10, -0.09, -0.92, -0.60, -1.20,
1971, -1.13, 0.24, -0.84, -0.24, 0.50, -1.57, 0.24, 1.55, 0.39, 0.58, -0.20, 0.60,
1972, 0.27, 0.32, 0.72, -0.22, 0.95, 0.88, 0.18, 1.32, -0.12, 1.09, 0.54, 0.19,
1973, 0.04, 0.85, 0.30, -0.54, -0.44, 0.39, 0.57, -0.06, -0.30, -1.24, -0.93, 0.32,
1974, 1.34, -0.14, -0.03, 0.51, -0.24, -0.14, -0.76, -0.64, 0.82, 0.49, -0.54, 1.50,
1975, 0.58, -0.62, -0.61, -1.60, -0.52, -0.84, 1.55, -0.26, 1.56, -0.54, 0.41, 0.00,
1976, -0.25, 0.93, 0.75, 0.26, 0.96, 0.80, -0.32, 1.92, -1.29, -0.08, 0.17, -1.60,
1977, -1.04, -0.49, -0.81, 0.65, -0.86, -0.57, -0.45, -0.28, 0.37, 0.52, -0.07, -1.00,
1978, 0.66, -2.20, 0.70, -1.17, 1.08, 1.38, -1.14, 0.64, 0.46, 1.93, 3.04, -1.57,
1979, -1.38, -0.67, 0.78, -1.71, -1.03, 1.60, 0.83, 0.96, 1.01, -0.30, 0.53, 1.00,
1980, -0.75, 0.05, -0.31, 1.29, -1.50, -0.37, -0.42, -2.24, 0.66, -1.77, -0.37, 0.78,
1981, 0.37, 0.92, -1.19, 0.36, 0.20, -0.45, 0.05, 0.39, -1.45, -1.35, -0.38, -0.02,
1982, -0.89, 1.15, 1.15, 0.10, -0.53, -1.63, 1.15, 0.26, 1.76, -0.74, 1.60, 1.78,
1983, 1.59, -0.53, 0.95, -0.85, -0.07, 0.99, 1.19, 1.61, -1.12, 0.65, -0.98, 0.29,
1984, 1.66, 0.72, -0.37, -0.28, 0.54, -0.42, -0.07, 1.15, 0.17, -0.07, -0.06, 0.00,
1985, -1.61, -0.49, 0.20, 0.32, -0.49, -0.80, 1.22, -0.48, -0.52, 0.90, -0.67, 0.22,
1986, 1.11, -1.00, 1.71, -0.59, 0.85, 1.22, 0.12, -1.09, -1.12, 1.55, 2.29, 0.99,
1987, -1.15, -0.73, 0.14, 2.00, 0.98, -1.82, 0.52, -0.83, -1.22, 0.14, 0.18, 0.32,
1988, 1.02, 0.76, -0.17, -1.17, 0.63, 0.88, -0.35, 0.04, -0.99, -1.08, -0.34, 0.61,
1989, 1.17, 2.00, 1.85, 0.28, 1.38, -0.27, 0.97, 0.01, 2.05, -0.03, 0.16, -1.15,
1990, 1.04, 1.41, 1.46, 2.00, -1.53, -0.02, 0.53, 0.97, 1.06, 0.23, -0.24, 0.22,
1991, 0.86, 1.04, -0.20, 0.29, 0.08, -0.82, -0.49, 1.23, 0.48, -0.19, 0.48, 0.46,
1992, -0.13, 1.07, 0.87, 1.86, 2.63, 0.20, 0.16, 0.85, -0.44, -1.76, 1.19, 0.47,
1993, 1.60, 0.50, 0.67, 0.97, -0.78, -0.59, -3.18, 0.12, -0.57, -0.71, 2.56, 1.56,
1994, 1.04, 0.46, 1.26, 1.14, -0.57, 1.52, 1.31, 0.38, -1.32, -0.97, 0.64, 2.02,
1995, 0.93, 1.14, 1.25, -0.85, -1.49, 0.13, -0.22, 0.69, 0.31, 0.19, -1.38, -1.67,
1996, -0.12, -0.07, -0.24, -0.17, -1.06, 0.56, 0.67, 1.02, -0.86, -0.33, -0.56, -1.41,
1997, -0.49, 1.70, 1.46, -1.02, -0.28, -1.47, 0.34, 0.83, 0.61, -1.70, -0.90, -0.96,
1998, 0.39, -0.11, 0.87, -0.68, -1.32, -2.72, -0.48, -0.02, -2.00, -0.29, -0.28, 0.87,
1999, 0.77, 0.29, 0.23, -0.95, 0.92, 1.12, -0.90, 0.39, 0.36, 0.20, 0.65, 1.61,
2000, 0.60, 1.70, 0.77, -0.03, 1.58, -0.03, -1.03, -0.29, -0.21, 0.92, -0.92, -0.58,
2001, 0.25, 0.45, -1.26, 0.00, -0.02, -0.20, -0.25, -0.07, -0.65, -0.24, 0.63, -0.83,
2002, 0.44, 1.10, 0.69, 1.18, -0.22, 0.38, 0.62, 0.38, -0.70, -2.28, -0.18, -0.94,
2003, 0.16, 0.62, 0.32, -0.18, 0.01, -0.07, 0.13, -0.07, 0.01, -1.26, 0.86, 0.64,
2004, -0.29, -0.14, 1.02, 1.15, 0.19, -0.89, 1.13, -0.48, 0.38, -1.10, 0.73, 1.21,
2005, 1.52, -0.06, -1.83, -0.30, -1.25, -0.05, -0.51, 0.37, 0.63, -0.98, -0.31, -0.44,
2006, 1.27, -0.51, -1.28, 1.24, -1.14, 0.84, 0.90, -1.73, -1.62, -2.24, 0.44, 1.34,
2007, 0.22, -0.47, 1.44, 0.17, 0.66, -1.31, -0.58, -0.14, 0.72, 0.45, 0.58, 0.34,
2008, 0.89, 0.73, 0.08, -1.07, -1.73, -1.39, -1.27, -1.16, 1.02, -0.04, -0.32, -0.28,
2009, -0.01, 0.06, 0.57, -0.20, 1.68, -1.21, -2.15, -0.19, 1.51, -1.03, -0.02, -1.93,
2010, -1.11, -1.98, -0.88, -0.72, -1.49, -0.82, -0.42, -1.22, -0.79, -0.93, -1.62, -1.85,
2011, -0.88, 0.70, 0.61, 2.48, -0.06, -1.28, -1.51, -1.35, 0.54, 0.39, 1.36, 2.52,
2012, 1.17, 0.42, 1.27, 0.47, -0.91, -2.53, -1.32, -0.98, -0.59, -2.06, -0.58, 0.17,
2013, 0.35, -0.45, -1.61, 0.69, 0.57, 0.52, 0.67, 0.97, 0.24, -1.28, 0.90, 0.95,
2014, 0.29, 1.34, 0.80, 0.31, -0.92, -0.97, 0.18, -1.68, 1.62, -1.27, 0.68, 1.86,
2015, 1.79, 1.32, 1.45, 0.73, 0.15, -0.07, -3.18, -0.76, -0.65, 0.44, 1.74, 2.24,
2016, 0.12, 1.58, 0.73, 0.38, -0.77, -0.43, -1.76, -1.65, 0.61, 0.41, -0.16, 0.48,
2017, 0.48, 1.00, 0.74, 1.73, -1.91, 0.05, 1.26, -1.10, -0.61, 0.19, -0.00, 0.88,
2018, 1.44, 1.58, -0.93, 1.24, 2.12, 1.09, 1.39, 1.97, 1.67, 0.93, -0.11, 0.61,
2019, 0.59, 0.29, 1.23, 0.47, -2.62, -1.09, -1.43, -1.17, -0.16, -1.41, 0.28, 1.20,
2020, 1.34, 1.26, 1.01, -1.02, -0.41, -0.15, -1.23, 0.12, 0.98, -0.65, 2.54, -0.30,
2021, -1.11, 0.14, 0.73, -1.43, -1.24, 0.77, 0.03, -0.28, -0.21, -2.29, -0.18, 0.29
)
#####
# Reformat and add a center date
month_centers <- data.frame(Month = month.abb,
Month_num = str_pad(seq(1,12,1), width = 2, pad = "0", side = "left"))
# Reformat
nao_long <- nao_idx %>%
pivot_longer(names_to = "Month", values_to = "nao_z", cols = Jan:Dec) %>%
left_join(month_centers, by = "Month") %>%
mutate(Date = as.Date(str_c(YEAR, "-", Month_num, "-15")),
ribbon_fill = ifelse(nao_z > 0,"above", "below"),
ribbon_above = ifelse(nao_z > 0, nao_z, NA),
above_ymin = ifelse(nao_z > 0, 0, NA),
ribbon_below = ifelse(nao_z < 0, nao_z, NA),
below_ymin = ifelse(nao_z < 0, 0, NA)
)
# Make figure
ggplot(filter(nao_long, YEAR >= 1982),
aes(x = Date, y = nao_z)) +
geom_ribbon(aes(ymin = above_ymin, ymax = ribbon_above, fill = "Positive Forcing"), alpha = 0.8) +
geom_ribbon(aes(ymin = below_ymin, ymax = ribbon_below, fill = "Negative Forcing"), alpha = 0.8) +
scale_fill_manual(values = c("Positive Forcing" = col_rdbu[1], "Negative Forcing" = col_rdbu[8])) +
geom_line(size = 0.25, color = "gray20") +
labs(x = "Date", y = "NAO", title = "NAO Index",
caption = "Years overlapping OISST data displayed", fill = "") +
theme_gmri() +
theme(legend.position = "bottom")
Looking at the NAO index, we found there was also not a strong relationship to our region’s Sea surface temperature anomalies. This lends strength to the argument that recent changes are independent of normal climate oscillations, but is not the end of the story.
Daily Temperature Record
The Gulf of Maine has a regular seasonal cycle hitting its lowest temperatures in March, and experiencing its peak highs in August. The difference between these two seasons is ~\(20-25^oF\) or ~\(11-12^oC\) (eyeballing it here).
hide
# Last Calendar Year
this_yr <- region_hw %>%
filter(year(time) == 2021)
# Set colors by name
color_vals <- c(
"Sea Surface Temperature" = "royalblue",
"Heatwave Event" = "darkred",
"MHW Threshold" = "coral3",
"Daily Climatology" = "gray30")
# Set the label with degree symbol
ylab <- "Sea Surface Temperature"
# Plot the last 365 days
hw_temp_p <- ggplot(this_yr, aes(x = time)) +
geom_segment(aes(x = time, xend = time, y = seas, yend = sst),
color = "royalblue", alpha = 0.25) +
geom_segment(aes(x = time, xend = time, y = mhw_thresh, yend = hwe),
color = "darkred", alpha = 0.25) +
geom_line(aes(y = sst, color = "Sea Surface Temperature")) +
geom_line(aes(y = hwe, color = "Heatwave Event")) +
geom_line(aes(y = mhw_thresh, color = "MHW Threshold"), lty = 3, size = .5) +
# geom_textpath(aes(y = mhw_thresh, color = "MHW Threshold"), label = "Heatwave Threshold", hjust = 0.8, lty = 3) +
#geom_line(aes(y = seas, color = "Daily Climatology"), lty = 2, size = .5) +
geom_textpath(aes(y = seas, color = "Daily Climatology"), label = "Climatology Mean", hjust = 0.5, lty = 2) +
scale_color_manual(values = color_vals) +
scale_x_date(date_labels = "%b", date_breaks = "1 month") +
scale_y_continuous(sec.axis = sec_axis(trans = ~as_fahrenheit(., data_type = "temperature"),
labels = number_format(suffix = " \u00b0F")),
labels = number_format(suffix = " \u00b0C")) +
theme(legend.title = element_blank(),
legend.position = "none") +
labs(x = "",
y = ylab)
# Plot the last 365 days - anomaly scale
linetype_key <- c(
"Sea Surface Temperature Anomaly" = 1,
"Heatwave Event" = 1,
"MHW Threshold" = 3,
"Daily Climatology" = 2)
# Same Plot as Anomalies
hw_anom_p <- this_yr %>%
mutate(
sst = sst,
seas = seas,
sst_anom = sst_anom,
mhw_thresh = mhw_thresh,
anom_thresh = mhw_thresh - seas,
anom_hwe = hwe - seas) %>%
ggplot(aes(x = time)) +
geom_segment(aes(x = time, xend = time,
y = 0, yend = sst_anom),
color = "royalblue", alpha = 0.25) +
geom_segment(aes(x = time, xend = time,
y = anom_thresh, yend = anom_hwe),
color = "darkred", alpha = 0.25) +
geom_line(aes(y = sst_anom, color = "Sea Surface Temperature Anomaly")) +
geom_line(aes(y = anom_hwe, color = "Heatwave Event")) +
geom_line(aes(y = anom_thresh, color = "MHW Threshold"), lty = 3, size = .5) +
geom_line(aes(y = 0, color = "Daily Climatology"), lty = 2, size = .5) +
# geom_textpath(aes(y = 0, color = "Daily Climatology"),
# label = "Climatology Mean", hjust = 0.5, lty = 2, show.legend = FALSE) +
scale_color_manual(values = color_vals) +
scale_linetype_manual(values = linetype_key, guide = "none") +
scale_x_date(date_labels = "%b", date_breaks = "1 month") +
guides(color = guide_legend(override.aes = list(linetype = linetype_key), nrow = 2)) +
theme(legend.title = element_blank(),
legend.position = "top") +
scale_y_continuous(sec.axis = sec_axis(trans = ~as_fahrenheit(., data_type = "anomalies"),
labels = number_format(suffix = " \u00b0F")),
labels = number_format(suffix = " \u00b0C")) +
labs(x = "",
y = "Temperature Anomaly",
caption = paste0("Climate reference period : 1982-2011"))
# Show figure
hw_temp_p / hw_anom_p
2021 Heatwave Events
hide
# Get Event numbers
events_21 <- this_yr %>% drop_na(mhw_event_no) %>% distinct(mhw_event_no) %>% pull(mhw_event_no)
events_num <- length(events_21)
# Pull their data
events_dat <- filter(region_hw, mhw_event_no %in% as.character(events_21),
yr <= 2021)
# Days in HW state
days_hw <- nrow(events_dat)
# How long on average
avg_duration <- round(days_hw / events_num, digits = 0)
During 2021 there were 5 main marine heatwave events, including one that carried over from the previous year and one that has continued into 2022. The average duration for these events was 72 days (excludes ongoing days in 2022).
Heatwave Seasonality
Over the last ~40 years the frequency, duration, and intensity of heatwave events has increased. This year was no exception with 358 days meeting the criteria for heatwave status.
hide
# Set new axis dimensions, y = year, x = day within year
# use a flate_date so that they don't stair step
base_date <- as.Date("2000-01-01")
grid_data <- region_hw %>%
mutate(year = year(time),
yday = yday(time),
flat_date = as.Date(yday-1, origin = base_date))
# Polar plot comparing when each year experienced peak heatwave conditions
polar_dat <- grid_data %>% filter(year %in% c("2012", "2021")) %>%
mutate(hw_line = ifelse(mhw_event == TRUE, sst_anom, NA),
hw_mag = ifelse(year == "2012", hw_line * -1, hw_line))
hide
# Things to highlight:
# when the heatwaves ocurred, how they overlapped
# severity?
# Side by side - need data separately
polar_dat %>%
ggplot(aes(x = flat_date, xend = flat_date,
y = sst_anom, yend = 0, color = mhw_event)) +
geom_segment(alpha = 0.9, show.legend = F) +
geom_textpath(aes(x = flat_date, y = mhw_thresh - seas),
label = "Heatwave Threshold", linetype = 1, color = "black", hjust = 0.85, straight = T) +
scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = expansion(mult = c(0,0))) +
scale_y_continuous(sec.axis = sec_axis(trans = ~as_fahrenheit(., data_type = "anomalies"),
labels = number_format(suffix = " \u00b0F")),
labels = number_format(suffix = " \u00b0C"),
expand = expansion(mult = c(0,0))) +
facet_wrap( ~ year, nrow = 2) +
scale_color_manual(values = c("TRUE" = "darkred", "FALSE" = "black")) +
# scale_color_gmri() +
# theme_gmri() +
labs(x = "", y = "Temperature Anomalies", fill = "Year")
Following Each Heatwave
Do some breakdowns of how each heatwave looked: - where was heat concentrated - how long did it last - was anything unique about it?
Heatwave Frequency Heatmap
hide
# Prep the legend title
guide_lab <- "Sea Surface Temperature Anomaly \u00b0C"
# Set palette limits to center it on 0 with scale_fill_distiller
limit <- max(abs(grid_data$sst_anom)) * c(-1, 1)
# Assemble heatmap plot
heatwave_heatmap <- ggplot(grid_data, aes(x = flat_date, y = year)) +
# background box fill for missing dates
geom_rect(xmin = base_date, xmax = base_date + 365,
ymin = min(grid_data$year) - .5, ymax = max(grid_data$year) + .5,
fill = "gray75", color = "transparent") +
# tile for sst colors
geom_tile(aes(fill = sst_anom)) +
# points for heatwave events
geom_point(data = filter(grid_data, mhw_event == TRUE),
aes(x = flat_date, y = year), size = .25) +
scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0)) +
scale_y_continuous(limits = c(1980.5, 2021.5), expand = c(0,0)) +
labs(x = "",
y = "",
"\nClimate reference period : 1982-2011") +
#scale_fill_gradient2(low = "blue", mid = "white", high = "red") +
scale_fill_distiller(palette = "RdYlBu", na.value = "transparent",
limit = limit) +
#5 inches is default rmarkdown height for barheight
guides("fill" = guide_colorbar(title = guide_lab,
title.position = "right",
title.hjust = 0.5,
barheight = unit(4.8, "inches"),
frame.colour = "black",
ticks.colour = "black")) +
theme_classic() +
theme(legend.title = element_text(angle = 90))
# Assemble pieces
heatwave_heatmap
Temperature Horizons
hide
# make slight modification to grid_data
hori_data <- grid_data %>%
mutate(year = factor(year),
year = fct_rev(year))
# cut out outliers and set cutpoints
cutpoints <- hori_data %>%
mutate(
outlier = between(
sst_anom,
quantile(sst_anom, 0.25, na.rm=T)-
1.5*IQR(sst_anom, na.rm=T),
quantile(sst_anom, 0.75, na.rm=T)+
1.5*IQR(sst_anom, na.rm=T))) %>%
filter(outlier)
# Set origin
ori <- 0
# Manually pick cut points
sca <- seq(-4, 4, length.out = 9)
sca_bins <- str_c(sca[1:(length(sca)-1)], " to ", sca[2:length(sca)], "\u00b0C")
sca_labels <- rev(sca_bins)
# Function to plot one or more years
anom_horizon_plot <- function(hori_data, origin = ori, scale_cutpoints = sca, labels = sca_labels){
ggplot(hori_data) +
geom_horizon(aes(flat_date,
sst_anom,
fill = ..Cutpoints..),
origin = ori,
horizonscale = sca) +
scale_fill_hcl(palette = 'RdBu', reverse = T, labels = sca_labels) +
facet_grid(year~.) +
theme(
panel.grid = element_blank(),
panel.spacing.y=unit(0, "lines"),
strip.text.y = element_text(size = 7, angle = 0, hjust = 0),
axis.text.y = element_blank(),
axis.title.y = element_blank(),
axis.ticks.y = element_blank(),
panel.border = element_blank(),
legend.position = "left"
) +
scale_x_date(expand=c(0,0),
date_breaks = "1 month",
date_labels = "%b") +
labs(x = '',
fill = "Temperature Anomaly \u00b0C")
}
# Plot Horizon Plot Using All Years
anom_horizon_plot(hori_data = filter(hori_data, yr <= 2021),
origin = ori,
scale_cutpoints = sca) +
labs(title = "Temperature Anomaly Horizons",
subtitle = "How deviation from the norm has become the new-normal",
caption = "Climatatological Reference Period: 1982-2011")
Comparing the Gulf’s Two Hottest Years
Early into 2021 it was apparent that the year was on-par with the previous title-holder for warmest year on record.
hide
# Reset the outliers
# cut out outliers and set cutpoints
cutpoints_2 <- hori_data %>%
filter(year %in% c("2012", "2021")) %>%
mutate(
outlier = between(
sst_anom,
quantile(sst_anom, 0.25, na.rm=T)-
1.5*IQR(sst_anom, na.rm=T),
quantile(sst_anom, 0.75, na.rm=T)+
1.5*IQR(sst_anom, na.rm=T))) %>%
filter(outlier)
# 2012
hori_2012 <- hori_data %>%
filter(year == "2012") %>%
anom_horizon_plot(scale_cutpoints = cutpoints_2) +
facet_wrap(~year, strip.position = "top")
# 2021
hori_2021 <- hori_data %>% filter(year == "2021") %>%
anom_horizon_plot(scale_cutpoints = cutpoints_2) +
facet_wrap(~year, strip.position = "top") +
theme(plot.caption = element_text(colour = "gray25", size = 6)) +
labs(caption = "Climatatological Reference Period: 1982-2011")
# Put them together as one figure
comparison_horizons <- hori_2012 / hori_2021 + plot_layout(guides = "collect")
# Do a bunch of formatting
comparison_horizons &
theme(
#Legend
legend.position = "top",
legend.key.height = unit(.5, "lines"),
legend.key.width = unit(3.75, "lines")) &
guides(fill = guide_legend(title = "Temperature Anomaly",
title.hjust = 0.5,
title.position = "top",
label.position = "bottom",
nrow = 1,
byrow = T, reverse = T)) &
plot_annotation(title = "Comparing the Gulf of Maine's Two Hottest Years",
subtitle = "Duration of temperature anomaly horizons, colored by their strength")
grid::grid.raster(lab_logo, x = 0.08, y = 0.04, just = c('left', 'bottom'), width = unit(0.8, 'inches'))
Streamplot
hide
# Set Factor Order for Bins
bin_order <- c(
"Less than -1",
"-1 to -0.5",
"Within 0.5",
"0.5 to 1",
"Greater than 1")
# Set Colors
bin_colors <- RColorBrewer::brewer.pal(n = length(bin_order),
name = "RdBu")
bin_colors <- rev(bin_colors)
# Count Days in Bins
day_types <- region_hw %>%
mutate(
`Less than -1` = if_else(sst_anom <= -1, TRUE, FALSE),
`-1 to -0.5` = if_else(between(sst_anom, -1, -0.5), TRUE, FALSE),
`Within 0.5` = if_else(between(sst_anom, -0.5, 0.5), TRUE, FALSE),
`0.5 to 1` = if_else(between(sst_anom, 0.5, 1), TRUE, FALSE),
`Greater than 1` = if_else(sst_anom >= 1, TRUE, FALSE)) %>%
group_by(yr) %>%
summarise(
`Less than -1` = sum(`Less than -1`, na.rm = T),
`-1 to -0.5` = sum(`-1 to -0.5`, na.rm = T),
`Within 0.5` = sum(`Within 0.5`, na.rm = T),
`0.5 to 1` = sum(`0.5 to 1`, na.rm = T),
`Greater than 1` = sum(`Greater than 1`, na.rm = T),
.groups = "drop") %>%
pivot_longer(cols = "Less than -1":"Greater than 1", names_to = "Anomaly Strength", values_to = "day_total", values_drop_na = FALSE) %>%
mutate(day_total = if_else(is.na(day_total), 0L, day_total),
`Anomaly Strength` = factor(`Anomaly Strength`, levels = bin_order))
# Plot the Streamplot
hw_streamplot <- ggplot(day_types, aes(yr, y = day_total, fill = `Anomaly Strength`)) +
geom_stream(type = "proportion") +
#geom_segment(y = .5, yend = 0.5, linetype = 3, color = "gray50", x = 1981, xend = 2021, size = 0.25) +
scale_x_continuous(breaks = seq(1980,2030, by = 10),
minor_breaks = seq(1975, 2030, by = 5),
expand = c(.03,0)) +
scale_y_continuous(labels = scales::percent_format(), expand = c(0,0)) +
scale_fill_manual(values = setNames(bin_colors, bin_order)) +
labs(y = "", x = "", fill = "",
title = "A Warmer Gulf of Maine",
subtitle = "Fractions of Each Year Experienced as Hot/Cold Anomalies",
caption = "1981 not a complete years in the data*") +
theme_minimal(base_size = 10) +
theme(
axis.line.x = element_line(),
axis.ticks.x = element_line(),
axis.line.y = element_line(),
axis.ticks.y = element_line(),
legend.position = "top",
legend.key.height = unit(.5, "lines"),
legend.key.width = unit(3.75, "lines"),
panel.grid.minor = element_blank(),
panel.grid = element_line(size = 0.2),
plot.margin = unit(c(.5,1,.3,.5), "cm"),
plot.title.position = "plot",
plot.title = element_text(face = "bold", size = 14),
plot.subtitle = element_text(size = 9),
plot.caption = element_text(size = 7.2, margin = margin(t = 20)),
legend.title = element_text(size = 9),
legend.text = element_text(size = 7.5),
legend.margin = margin(c(7,0,0,0)),
legend.justification = "center") +
guides(fill = guide_legend(title = "Anomaly Strength \u00b0C",
title.hjust = 0.5,
title.position = "left",
label.position = "bottom",
nrow = 1,
byrow = T))
# Show figure
hw_streamplot
grid::grid.raster(lab_logo, x = 0.08, y = 0.04, just = c('left', 'bottom'), width = unit(0.8, 'inches'))
Characteristics of the Hottest Years
Multivariate figure comparing characteristics of both years: 2012 & 2021
hide
# Use all years so they can be ranked independently
# Metrics?
year_mets <- hori_data %>%
filter(year != "2022") %>%
group_by(year) %>%
summarise(
total_days = n(),
`Average Temperature` = round(mean(sst, na.rm = T), 2),
`Temperature Above Average` = round(mean(sst_anom, na.rm = T), 2),
`Cumulative Degrees Above Average` = round(sum(sst_anom, na.rm = T), 0),
hw_days = sum(mhw_event, na.rm = T),
hw_events = n_distinct(mhw_event_no),
peak_temp = max(sst),
peak_anom = max(sst_anom)) %>%
mutate(
hw_day_pct = round((hw_days/total_days) * 100, 2),
`Average Heatwave Length` = round(hw_days/hw_events, 1)
) %>%
pivot_longer(names_to = "Var", values_to = "Metric", cols = `Average Temperature`:`Average Heatwave Length`)
# Snipe out the top5 for each
year_met_ranks <- year_mets %>%
filter(Var %in% c("Average Temperature", "Temperature Above Average", "Average Heatwave Length", "Cumulative Degrees Above Average")) %>%
group_by(Var) %>%
slice_max(n = 5, order_by = Metric) %>%
ungroup() %>%
mutate(
yr_focus = ifelse(year == "2021", TRUE, FALSE),
year = reorder_within(year, Metric, Var))
year_met_ranks %>%
ggplot(aes(year, Metric, color = yr_focus)) +
geom_col(aes(fill = yr_focus), show.legend = F) +
geom_label(aes(label = Metric), show.legend = FALSE) +
facet_wrap(~Var, scales = "free", ncol = 2) +
coord_flip() +
scale_x_reordered() +
scale_y_continuous(expand = expansion(mult = c(0, .2))) +
labs(y = "",
x = NULL,
title = "How 2021 Stacks Up:",
subtitle = "Ranking Characteristics of Acute and Prolonged Warming Events") +
theme(plot.title = element_text(hjust = 0)) +
scale_color_gmri() +
scale_fill_gmri() +
theme(panel.border = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.major. = element_line(linetype = 3, color = "gray50"))
Maps
- Where was 2021 the 1st, 2nd, third hottest year
hide
# Load the rankings done in python, do some checking
ranks_path <- str_c(oisst_path, 'temp_rankings/yrly_ranks_1982to2021.nc')
# Pull this year
ranks_21 <- stack(ranks_path)[["X2021"]]
# Reclassify Ranks
reclass_df <- c(
0, 1, -1, # Coldest Year
1, 2, -2, # Second Coldest
2, 3, -3, # Third Coldest
3, 37, NA, # Things to Mask
37, 38, 3, # Third Warmest
38, 39, 2, # Second Warmest
39, 40, 1 # Warmest
)
# make reclassification matrix
reclass_m <- matrix(reclass_df, ncol = 3, byrow = TRUE)
# reclassify - masking middle values handled with reclassification
ranks_reclass <- reclassify(ranks_21, reclass_m)
# Make Stars for ggplot
ranks_st <- st_as_stars(rotate(ranks_reclass))
# Mutate to get levels as factors for legend
ranks_st <- ranks_st %>%
mutate(
temp_rank = case_when(
X2021 == -1 ~ "Coldest",
X2021 == -2 ~ "2nd Coldest",
X2021 == -3 ~ "3rd Coldest",
X2021 == 3 ~ "3rd Warmest",
X2021 == 2 ~ "2nd Warmest",
X2021 == 1 ~ "Warmest",
TRUE ~ "Not a Record Year"),
temp_rank = fct_drop(temp_rank),
temp_rank = factor(temp_rank, levels = c(
"Warmest", "2nd Warmest", "3rd Warmest", "Not a Record Year",
"3rd Coldest", "2nd Coldest", "Coldest"))) %>%
select(-X2021)
# Make Map
ggplot() +
geom_stars(data = ranks_st) +
geom_sf(data = world_sf, fill = "gray90", size = .25, alpha = 0.2) +
scale_fill_brewer(palette = "RdBu", na.value = "transparent") +
scale_x_continuous(breaks = seq(-180, 180, 90),
expand = expansion(mult = c(0,0))) +
scale_y_continuous(breaks = seq(-90, 90, 30),
expand = expansion(mult = c(0,0))) +
labs(title = "Record Setting Temperatures of 2021",
subtitle = "Regions of the planet that experienced record setting temperatures",
fill = "",
caption = "Temperature rankings cover period of 1982-2021") +
map_theme
- Where is it warming fastest
hide
# Load the LME timeseries to calculate rates of warming
lme_names <- get_region_names("lme") # names
lme_paths <- get_timeseries_paths("lme", mac_os = "mojave") # paths
# Data
lme_oisst <- map(lme_names, ~ read_csv(lme_paths[[.x]][["timeseries_path"]], col_types = cols(), guess_max = 1e6))
lme_oisst <- setNames(lme_oisst, lme_names)
# Add Gulf of Maine Timeseries to the List with LME's
lme_oisst[["gulf_of_maine"]] <- region_timeseries
# Make sure dates are formatted
lme_oisst <- map(lme_oisst, ~ mutate(.x, time = as.Date(time)))
#### Warming Rates ####
lme_8221 <- map(lme_oisst, ~ mutate(.x, year = year(time)) %>% filter(year %in% c(1982:2021)))
# Get warming rates
lme_rates <- map_dfr(lme_8221, function(lme_data){
# Uses area weighted SST
lme_yearly <- group_by(lme_data, year) %>%
summarise(sst = mean(area_wtd_sst, na.rm = T))
yrly_warming <- lm(sst ~ year, data = lme_yearly)
mod_details <- broom::tidy(yrly_warming)
yrly_rate <- mod_details[2,"estimate"]
names(yrly_rate) <- "annual_warming_rate_C"
yrly_rate
}, .id = "Region") %>%
arrange(desc(annual_warming_rate_C))
# Pull out the top 13
rates_ordered <- lme_rates %>%
mutate(`Warming Rank` = row_number(),
`Warming Rate F` = as_fahrenheit(annual_warming_rate_C, data_type = "anomalies"),
Region = str_replace_all(Region, "_", " "),
Region = str_to_title(Region)) %>%
select(Region, `Warming Rank`, `Warming Rate C` = annual_warming_rate_C, `Warming Rate F`)
# symbol for degree: \u00b0
# Make table of top 10*
rates_ordered %>%
slice(1:10) %>%
mutate_if(is_numeric, round, 3) %>%
gt(rowname_col = "Region") %>%
tab_stubhead(label = "Region") %>%
tab_header(
title = md("**Warming Rate Rankings of Large Marine Ecosystems**")) %>%
tab_source_note(source_note = md("*Data Source: NOAA OISSTv2 Daily Sea Surface Temperature Data. Temperatures adjusted for latitudinal distotion.*") ) %>%
tab_source_note(source_note = md("*Gulf of Maine not a true LME*") )
| Warming Rate Rankings of Large Marine Ecosystems | |||
|---|---|---|---|
| Region | Warming Rank | Warming Rate C | Warming Rate F |
| Baltic Sea | 1 | 0.047 | 0.085 |
| Gulf Of Maine | 2 | 0.044 | 0.080 |
| Black Sea | 3 | 0.044 | 0.079 |
| Scotian Shelf | 4 | 0.041 | 0.073 |
| Iceland Shelf And Sea | 5 | 0.038 | 0.069 |
| Northeast Us Continental Shelf | 6 | 0.038 | 0.068 |
| North Sea | 7 | 0.035 | 0.063 |
| Norwegian Sea | 8 | 0.032 | 0.057 |
| Sea Of Japan | 9 | 0.031 | 0.056 |
| Mediterranean Sea | 10 | 0.030 | 0.055 |
| Data Source: NOAA OISSTv2 Daily Sea Surface Temperature Data. Temperatures adjusted for latitudinal distotion. | |||
| Gulf of Maine not a true LME | |||
hide
#
rates_path <- str_c("{box_root}RES_Data/OISST/oisst_mainstays/warming_rates/annual_warming_rates{trend_start_year}to{trend_end_year}.nc")
# Get paths to the LME's and their timeseries
lme_paths <- gmRi::get_timeseries_paths(region_group = "lme", mac_os = "mojave")
- Comparing Regions in Table
Biological Events
- basking sharks
- hammerheads
- blue crab
A work by Adam A. Kemberling
Akemberling@gmri.org